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Quantum Monte Carlo results for the specific heat c of the two dimensional Hubbard model are 
presented. At half-filling it was observed that c ~ at very low temperatures. Two distinct 
features were also identified: a low temperature peak related to the spin degrees of freedom and a 
higher temperature broad peak related to the charge degrees of freedom. Away from half-filling the 
spin induced feature slowly disappears as a function of hole doping while the charge feature moves 
to lower temperature. A comparison with experimental results for the high temperature cuprates is 
discussed. 

PACS numbers; 65.40.-|-g, 65.50.-fm, 75.40.Cx, 75.40.Mg 



I. INTRODUCTION 



II. HALF-FILLING 



The Hubbard model is among the simplest Hamilto- 
nians that describe the behavior of correlated electrons. 
Specially since the discovery of high temperature super- 
conducting materials, considerable attention has been de- 
voted to this model and significant progress was achieved 
in understanding its ground state properties, particularly 
at half-filling, although superconductivity is still elusive 
0. Static and dynamical spin correlations, the opti- 
cal conductivity and other observables have been stud- 
ied in detail, [y However, not much attention has been 
devoted to its thermodynamical properties despite the 
large amount of experimental specific heat measurements 
available for the cuprates. The aim of this paper is to fill 
that void and to present a systematic study of the spe- 
cific heat of the two dimensional (2D) Hubbard model for 
different couplings U/t, dopings and temperatures. To 
achieve that goal Quantum Monte Carlo (QMC) tech- 
niques are used. 

The Hubbard Hamiltonian is given by 



h.c.) 



(1) 



where creates an electron at site i with spin projec- 
tion (T, riio- is the number operator, the sum (ij) runs 
over pairs of nearest neighbor lattice sites, U is the on- 
site Coulombic repulsion, t the nearest neighbor hopping 
amplitude, and /i the chemical potential. In the following 
t ~ \ will be used as the unit of energy. The boundary 
conditions are periodic. 



The computational calculation of the specific heat c is 
not simple. In principle, c is given by the derivative of 
the energy E (defined as E = {H)/N, with N being the 
number of sites) with respect to the temperature T at 
constant density. However, note that in dctcrminantal 
QMC simulations, which are set up in the grand canon- 
ical ensemble, the energy is a function of the chemical 
potential that has to be adjusted to keep the density (n) 
constant as the temperature changes. In other words, 
dE/dT must be calculated along lines of constant (n) in 
the T — II plane. In this framework the calculation of 
c cannot proceed using c ~ {H'^) — {H)'^, as when the 
number of particles is fixed. Another detail that is im- 
portant is the finite discretization of the derivatives along 
the lines of constant density. Naively, the ratio AE/AT, 
with AT very small, should be calculated. However, us- 




FIG. 1. Monte Carlo results for the energy £ on a 6 x 6 
cluster at half-filling and U = 8 (open circles). The low tem- 
perature polynomial fit is indicated by the dashed line, while 
the short-long dashed line indicates the high temperature fit. 
The solid line denotes the specific heat c. The low tempera- 
ture data that produce the spin peak are shown in the inset. 
The error bars are smaller than the size of the dots. 
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FIG. 2. Energy i5 as a function of temperature T at 
half-filling on 4 x 4, 6 x 6 and 8x8 clusters for a) [/ = 0, b) 
f/ = 4, c) [/ = 8 and d) f7 = 12. The error bars are smaller 
than the size of the dots. 

ing such a procedure the small statistical error in E intro- 
duces large errors in c. For that reason we have decided 
to calculate £'(T) (at fixed (n)) numerically as accurately 
as possible, and then fit the Monte Carlo points with a 
polynomial that smears out the small fluctuations in E. 
c is obtained by taking derivatives from this polynomial 
analytically. Motivated by the shape of the E vsT curve, 
different polynomials were used for the high and low tem- 
perature regimes. In Fig.l, the raw Monte Carlo data for 
E as a function of temperature corresponding to [/ = 8 
at half-filling on a 6 x 6 cluster are presented. Each data 
point was obtained by performing around 10,000 mea- 
surement sweeps. The dashed line indicates the low tem- 
perature fit by a polynomial of order 6 in T, while the 
short-long dashed line indicates the high temperature fit, 
in this case to a polynomial of order 4. T* is the temp- 

6x6, <n>=1.0 




FIG. 3. Specific heat c vs T at half-filling for different val- 
ues of U ranging from 2 to 12. The vertical axis for each 
coupling is shifted for clarity. 
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FIG. 4. a) Temperature Tspin peak, where the spin peak is 
located, as a function of U at half-filling. The dashed line indi- 
cates T — 2 J/3 (asymptotic result in the Heisenberg limit); b) 
Temperature Tcharge peak, where the charge peak is located, 
as a function of U. The dashed line indicates T = 0.24C/. 

erature where the two fits meet. Its value depends on 
the parameters U and (n) and is typically of the order 
of 1. In order to make a smooth connection of the two 
fits we included points below (above) T* for the high 
(low) temperature fit within a window ^ 0.2 centered at 
T*. The specific heat was obtained through the analytic 
derivative of the fitting polynomials, and it is also shown 
in Fig.l with a continuous line. The inset of the figure 
shows with more detail the low energy data that generate 
the low temperature peak in c (to be discussed later). 

An important issue in QMC simulations are finite size 
effects (FSE). Upon studying 4 x 4, 6 x 6 and 8x8 clus- 
ters, it was observed that the FSE in £' vs T are strong 
at very weak coupling but become negligible for U = 8 
or larger. In Fig. 2, the energy of the different clusters 
for U — 0, 4:, 8 and 12 is shown. Since the FSE are 
small we decided that results on 6 x 6 clusters are repre- 
sentative of the physical behavior analyzed in this study 
and, thus, this is the lattice size that we have used in the 
remaining of the paper. In Fig. 3, c vs T at half-filling 
for different values of U is shown. There are two impor- 
tant features in these curves: 1) A low temperature peak 
that appears when the low lying spin states are excited, 
and 2) a higher temperature peak which appears when 
states in the upper Hubbard band are excited. In the 
weak coupling regime the low temperature peak moves 
to slightly higher temperature as U increases, reaching 
a turning point at C/ « 7 where the peak is at T = 0.3. 
For U > 7 the peak slowly moves to lower temperatures, 
as U grows. This indicates the beginning of the strong 
coupling regime since it is well known that for large val- 
ues of U the Hubbard and the t — J models have similar 
behaviors and the coupling constants are related through 
J = At'^/U. Numerical studies on the t — J model have 
indicated that at half-filling (Heisenberg limit) the peak 
in c appears at T « 2 J/3 ||] which in terms of U corre- 
sponds to T « 8i^/3[/. Thus, when this regime is reached 
we expect the peak to move to lower temperature with 
increasing U. The position of the peak as a function of 
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FIG. 5. a) Energy vs at half- filling for different values 
of U\ b) Coefficient 5 of the low temperature fit c ~ 5T^ , 
as a function of U at half-filling. Numerical and analytical 
values for the Ifeisenberg limit {U — oo) are indicated. The 
circle corresponds to the mean field result of Ref. [6] , and the 
diamond denotes the numerical result of Ref. [8] . 

U/t is shown in Fig. 4. a, where the dashed line indicates 
T = 2J/3. The asymptotic behavior is reached for U > 
10. 

The broad high temperature peak moves to higher tem- 
perature as U increases as expected since its presence cor- 
responds to the excitation of states across the gap that 
grows with U. In Fig.4.b the position of this peak is 
shown as a function of U. For U > 7 the dependence of 
the position of the peak with U becomes approximately 
linear, and it is given by 0.24C/. A spin-density-wave 
mean field calculation of the gap A as a function oiU, at 
large U, gives the result A ~ 0.48J7. Apparently, quan- 
tum fluctuations reduce the size of the gap. Note that 
in Fig. 3 it can be observed that the minimum in c be- 
tween the two peaks becomes deeper as U increases and 
the charge peak increases its width. 

In previous work the specific heat for the half-filled 
Hubbard model in one dimension has been evaluated. 

We found that the qualitative behavior in two and 
one dimensions is similar regarding the existence and cou- 
pling dependence of the two peaks. However, the follow- 
ing differences were observed: 1) According to Ref. ||] 
the two peaks can be resolved for U > A while here we 
were able to identify the two peaks already at U = 2. 
The fact that only one maximum is observed in Ref. |^] 
in the strong coupling regime is due to the small T in- 
terval considered in their study; 2) According to Ref. 
the maximum in c associated with the spin excitations 
moves to lower temperatures as U increases in weak cou- 
pling while in our 2D study the opposite behavior was 
found. 

Another important feature observed here at half-filling 
is that at low temperatures the specific heat follows 
c = (5r^, i.e., the behavior predicted by spin- wave cal- 
culations. 1^ In Fig. 5. a we show the energy as a function 
of for different values of U showing that linear behav- 
ior occurs for T < 0.3. The value of S depends on U, and 
it decreases as the coupling increases. For large U the 
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FIG. 6. c vs T for i7 = 2 to 5 (continuous lines) and 
for ?7 = 6 to 12 (dashed lines). All the curves intersect at 
T-i ~ 1.6, while those corresponding to weak coupling (con- 
tinuous lines) also intersect at T\ ~ 0.6. 

limiting value S ~ 1.1 is obtained in good agreement 
with the reported value for the Heisenberg model. 
A slave boson mean field theory (SBMFT) calculation 
provided a value of 5 = 1.3 ± 0.05 |^ while a numerical 
study obtained 5 = 1.1 ± 0.2. g The behavior of S vs U 
is shown in Fig.5.b. 

In Fig. 6, c vs T for several values of U ranging from 2 
to 12 are presented. These are the same curves that were 
shown in Fig. 3 but now using common vertical units. It 
is interesting to observe that all the curves intersect at 
r=1.6±0.2. If only small values of the coupling are 
considered, i.e. U ranging from 2 to 5, the curves cross 
also at Ti = 0.6 in addition to T2 = 1.6. This behavior 
was predicted by VoUhardt |^ and was observed in the 
paramagnetic phase of the infinite dimensional Hubbard 
model for 0<U < 2.5. 



III. FINITE HOLE DENSITY 

To compare our results with those of the superconduct- 
ing cuprates it is important to study the specific heat as 
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FIG. 7. Energy £ as a function of temperature T on 4 x 4 
and 6x6 clusters and density {n) = 0.75 for &) U = Q, b) 
)7 = 4, c) (7 = 8 and d) [7 = 12. 
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FIG. 8. c vs T at density a) (n> = 0.75 and b) (n) = 0.50 
for [/ = 0, 4, 8 and 12, on a 6 x 6 cluster. 

a function of hole doping. As remarked before, many 
experimental measurements of the specific heat for high 
temperature cuprates are available. In general it is very 
difficult to separate the electronic contribution to the spe- 
cific heat in the normal state from the phononic part. 
Also many experiments have been performed in the su- 
perconducting phase, where the existence of an intrinsic 
linear contribution to the specific heat would indicate 
the absence of a gap and thus non-conventional behav- 
ior. [|l^,^^ Since the superconducting phase cannot 
be reached in QMC simulations, our results will be com- 
pared with experiments performed in the normal state. 
For La2-xSrxCu04 it was observed in Ref. ||l^ that the 
linear term 7 of the specific heat in the normal phase in- 
creases with doping between x = 0.12 and 0.25. However, 
studies of the same material performed later showed 
that 7 increases with x for a; > 0.1 reaching a maximum 
value at optimal doping x « 0.15 and then decreasing in 
the overdoped regime. This behavior is in agreement with 
the Van Hove scenario where the density of states 
reaches a maximum at optimal doping. The behavior 
of 7 for a metal-insulator transition was also studied for 
Sri-xLttxTiO^ in Ref. |l^. They observed that 7 in- 
creases as the transition is approached from the metallic 
side. Through the relation 7 = m*jo/m, where 70 and m 
are the linear coefficient and the mass for free electrons, 
it was found that the effective mass of the quasiparticles 
m* increases as the transition is approached. Loram et 
al. ||l^ studied 7 as a function of doping at T = 280K 
in Y Ba2C'u30e-f-x- 7 appears to increase with doping 
reaching a plateau for x « 0.45. 

The first step to study numerically the specific heat at 
finite density is to analyze the finite size effects. They 
are stronger than at half-filling, but still moderate as can 
be seen in Fig. 7 where E vs T for U = 0, 4, 8 and 12 at 
(n) = 0.75 on 4 X 4 and 6x6 clusters is shown. Away from 
half-filling it was very difficult to obtain accurate results 
on 8 X 8 clusters at low temperature due to the well- 
known sign problem. However, since FSE are stronger 
in weak coupling and we have observed that for U — 0, 
where results can be obtained exactly, there is only a 
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FIG. 9. c vs T for f/ = 8 (solid line) on a 6 x 6 cluster at 
different densities. The dashed line indicates results for (7 = 
obtained using a 200 x 200 cluster. 

small difference between the 6x6 and 8x8 results, then 
as before, 6x6 lattices were used in our studies away 
from half-filling. 

In Fig. 8. a the specific heat as a function of T at 
(n) = 0.75 for different values of U is presented. It can be 
observed that the spin peak is substantially reduced com- 
pared with the results at half-filling, but it is still present 
in strong coupling for ?7 = 8 and 12 indicating the ex- 
istence of short range antiferromagnetic correlations. In 
weak coupling, i.e. for U = A, the spin feature has dis- 
appeared, and the curve is similar to the non-interacting 
one. The specific heat increases in the region where the 
minimum between the two peaks existed at half-filling. 
At quarter filling (Fig.S.b), c has a behavior that resem- 
bles free electrons independently of the value of U. Thus, 
here the electrons are approximately weakly interacting 
at all couplings. 

Let us consider in more detail the special case of U=8. 
This value of the coupling was selected since according 
to calculations of the optical conductivity it is suitable 
to reproduce some normal state experimental results [0. 
In Fig. 9 the specific heat as a function of temperature 
is presented for different values of the density (n). The 
continuous line indicates the results for (7 = 8 on a 6 x 6 
cluster while the dashed line denotes the non-interacting 
U ~ results on a 200 x 200 lattice. Such a large cluster 
in the non-interacting case was used to avoid finite size 
effects which are strong in this limit at the low temper- 
atures where the linear behavior occurs. Again it should 
be remarked that this problem occurs in weak coupling 
at very low T and, thus, our U = 8 results are not ex- 
pected to be contaminated by size effects. In Fig. 9 it can 
be seen that for J7 = 8 the intensity of the spin peak de- 
creases smoothly with doping. At 10% hole doping (i.e. 
(n) — 0.90) its intensity diminishes by 40%, a result in 
agreement with Ref. ^ where the t — J model was stud- 
ied. Note that for (n) ^ 0.8 the specific heat is almost 
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FIG. 10. c/T vs T for ?7 = 8 (solid line) on a 6 x 6 cluster at 
different densities. The dashed lines denote results for U — 
obtained using a 200 x 200 cluster. 

flat in a broad range of temperatures. Here it is difflcult 
to resolve the spin and charge peaks from the data. We 
expect that at this density or lower the spin correlations 
are no longer important, even those of short-range, in 
agreement with previous spectral function studies per- 
formed in the Hubbard model. |l^ Reducing further the 
density from (n) =0.75 to 0.5 a single peak structure that 
resembles the non-interacting specific heat curve becomes 
dominant. 

An important issue in this context is the calculation of 
couplings and densities where the system changes from 
insulator to metal. Metallic behavior is characterized in 
the specific heat by the existence of a linear coefficient 7. 
In two dimensions it was found that [[191 



7T - 



(2) 



with T2D positive in strong coupling. 

The experimentalists often present plots of c/T vs 
when addressing 7. Analogously, in Fig. 10 the continu- 
ous line denotes c/T vs T for U — 8 at different densities, 
while the dashed line indicates the non-interacting case. 
The lowest temperature that was confidently reached in 
this study away from half-filling is T = 0.3. It is clear 
that this temperature is too high to observe the linear be- 
havior in c/T since, according to Eq.(2), the slope of the 
curve has to be positive at very low temperature. Clearly, 
if the system behaves as a Fermi liquid a maximum has to 
appear in the curve at a lower temperature than reached 
in this study. The non-interacting results show indeed 
the linear behavior at very low temperatures. However, 
note that the value of 7 for non-interacting electrons is 
not much different from the value of c/T at the maxi- 
mum in Fig. 10 at all densities. Thus by extrapolating 
the U = 8 curves to zero we expect to obtain a good 
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FIG. 11. c/T vs (n) for U = 8 at different temperatures. 

approximation to the value of 7. However, since we can- 
not reach lower temperatures the existence of anomalous 
non-Fermi liquid behavior can certainly not be ruled out 
as remarked in Ref . ||] . 

In Fig. 11, c/T as a function of doping is presented at 
different temperatures. Notice that the lowest tempera- 
ture T = 0.5 shown in Fig. 11 corresponds to 2000i4r, if 
t = OAeV is used. This is much higher than T = 280K 
which is the highest temperature used in experiments. 
However, for T = 0.5 it was here observed that c/T in- 
creases with doping for (n) < 0.8 in agreement with some 
experimental results p7| , and the same behavior is ob- 
served at r = 1. For higher T the ratio c/T increases for 
increasing density (n) . 



IV. SUMMARY 

The specific heat of the two dimensional Hubbard 
model has been calculated for different couplings and 
electronic densities as a function of temperature. At 
half-filling and as the coupling U increases a low temper- 
ature peak associated with the spin degrees of freedom 
moves to lower temperatures, while a high temperature 
feature associated with the charge degrees of freedom 
moves to higher temperatures. At very low temperatures 
c « ST^ as predicted by spin-wave theory and S tends 
to the Heisenberg value {S « 1.1) for large coupling U. 
Away from half-filling we observed that the spin feature 
weakens with doping, and it disappears for (n) < 0.75 
working at C/ = 8. This suggests the absence of impor- 
tant antiferromagnetic correlations below that density. 
We were not able to reach temperatures low enough to 
decide whether the system is metallic or has anomalous 
behavior away from half-filling. However, by evaluating 
c/T we were able to make comparisons with experimen- 
tal results. At the lowest temperatures that we could 
reach we found that c/T increases with hole doping for 
(n) < 0.9. This behavior is similar to experimental re- 
sults for YBa2Cu30e+x- 
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